function eff = best_ans_cut(in)

mir_wide=in(1);%镜宽
mir_length=in(2);%镜长
mir_hight=in(3);%镜高

%% 已知参数
month=3;
D=daysact(datetime(2023,3,21),datetime(2023,month,21));
D=D+365*(D<0);%日期
t=9:1.5:15;%时刻
lat=39.4/180*pi; %纬度
H=3; %海拔 km
delta=asin(sin(2*pi*D/365)*sin(2*pi*23.45/360));%赤纬角
tower=[0,0,80];%吸收塔坐标

%镜子离散化
dx=20;dy=20;
[xx,yy]=meshgrid(-mir_wide/2:mir_wide/dx:mir_wide/2,-mir_length/2:mir_length/dy:mir_length/2);
[X,Y]=distribution(mir_wide);%计算各个镜子的坐标位置
mir_center=[X,Y,mir_hight*ones(size(X))];%计算镜子中心坐标

dHR=tower-mir_center;
V_out=normr(dHR);%光的反射向量
dHR=sqrt(sum(dHR.^2,2));%计算镜子到吸收塔的直线距离

nsigma=10;
ntau=10;
[sigma,tau]=meshgrid(linspace(4.65e-3/nsigma,4.65e-3,nsigma),linspace(0,2*pi,ntau)); % 在一个圆锥内产生均匀的100条光线组成光束
sigma=sigma(:);tau=tau(:);
V_light=[sin(sigma).*cos(tau),sin(sigma).*sin(tau),cos(sigma)]; % 光线关于地面坐标系下的方向矢量

eff_t=zeros(length(X),5); % 每一列记录第t时刻下每片镜子的光学效率
power_t=zeros(5,1); % 记录每个t时刻下整个系统的功率
dni_t=zeros(5,1); % 5个不同时刻的DNI值
%% 计算
for temp_t=1:length(t)
    ST=t(temp_t);
    omega=pi/12*(ST-12);
    alpha_s=asin(cos(delta)*cos(lat)*cos(omega)+sin(delta)*sin(lat)); % 太阳高度角   
    gamma_s=real(acos((sin(delta)-sin(alpha_s)*sin(lat))/cos(alpha_s)/cos(lat))); % 太阳方位角  
    gamma_s=(2*pi-gamma_s)*(ST>12)+gamma_s*(ST<=12);
    
    V_in=[cos(alpha_s)*sin(gamma_s),cos(alpha_s)*cos(gamma_s),sin(alpha_s)];%光的入射向量
    light_X=acos(V_in(3)); % 主轴光线与Z轴偏离角度，就是说要先绕X正向转这些角度才行
    light_Z=angle(V_in(1)+1i*V_in(2))+pi/2;
    martix_trans_light=martix_trans(light_X,light_Z); % 给出光线的旋转变换矩阵
    V_light_1=(martix_trans_light*V_light')';
    
    V_mid=normr(V_in+V_out);
    angle_X=acos(V_mid(:,3));
    angle_Z=angle(V_mid(:,1)+1i*V_mid(:,2))+pi/2;
    
    
    %各个效率计算
    eff_sb=zeros(size(X));
    eff_cos=V_mid*V_in';
    eff_at=0.99321-0.0001176*dHR+1.97e-8*dHR.^2;
    eff_trun=eff_sb;
    eff_ref=0.92;
    
    dni=1.366*(0.4237-0.00821*(6-H)^2+(0.5055+0.00595*(6.5-H)^2)*exp(-(0.2711+0.01858*(2.5-H)^2)/sin(alpha_s)));
    
    for i=1:length(X)
        H1_i=[xx(:),yy(:),zeros(length(xx(:)),1)];
        matrix_change_i=martix_trans(angle_X(i),angle_Z(i));
        H2_i=(matrix_change_i*H1_i')'+mir_center(i,:);
        difference=([X,Y]-[X(i),Y(i)]);
        distan=abs(difference(:,1)+1i*difference(:,2));%记录镜子之间的相互距离
        distan_sort=sort(distan); % 从小到大排序
        close_6=find(distan<=distan_sort(6) & distan>0);%选择距离最近的6个镜子
        shadow_index=[];%用于记录阴影遮挡镜子序号
        %% 寻找遮挡镜子
        for j=1:length(close_6) 
            norm_j=V_mid(close_6(j),:); % 镜的法矢量
            jiao_s=(dot(norm_j,mir_center(close_6(j),:))-H2_i*norm_j')/dot(norm_j,V_in); % 这是i镜上各个离散点沿着光射入的方向去找与j镜所在平面的交点参数s
            
            node_loc=H2_i+jiao_s*V_in; % 交点的全局坐标    
            trans_matrix_j=martix_trans(angle_X(close_6(j)),angle_Z(close_6(j))); %  % 给出j镜的旋转变换矩阵，
            node_in_j=(trans_matrix_j\(node_loc-mir_center(close_6(j),:))')'; % 求逆,给出交点在j镜中的局部坐标
            shadow_index=[shadow_index;setdiff(find(abs(node_in_j(:,1))<mir_wide/2 & abs(node_in_j(:,2))<mir_length/2 & jiao_s>0),shadow_index)]; % 判断交点是否在j内,不重复记录已经判定为遮挡交点 
        end
        num_shadow=length(shadow_index); % i镜上的离散点总的被遮挡个数,计算遮挡效率
        eff_sb(i)=1-num_shadow/(dx+1)/(dy+1);%更新阴影遮挡效率
        
        %%  截断效率
        V_light_2=[V_light_1;V_in]; % 把中心光线也加上，一共101条，下面看这些光线在各个离散点反射后有多少能够到达吸收器
        V_light_out=2*V_light_2*V_mid(i,:)'.*V_mid(i,:)-V_light_2; % 计算光束里面所有光线的反射方向
        n_jie=0;
        
        %计算光在第k离散点反射后光线是否与接收器圆柱相交
        for k=1:length(H2_i) % 每个离散点被光束照射
            if sum(k==shadow_index)==0 % 如果不是被遮挡点，做计算
                point_k=H2_i(k,:); % 当前反射点坐标
                s_light_2=(-(2*point_k(1)*V_light_out(:,1)+2*point_k(2)*V_light_out(:,2))-sqrt((2*point_k(1)*V_light_out(:,1)+2*point_k(2)*V_light_out(:,2)).^2-....
                4*(V_light_out(:,1).^2+V_light_out(:,2).^2)*(point_k(1)^2+point_k(2)^2-3.5^2)))./(2*(V_light_out(:,1).^2+V_light_out(:,2).^2));
                
            
                index_real=real(s_light_2)==s_light_2; % 哪些光线交点有解
                z_jiao2=point_k(3)+s_light_2(index_real).*V_light_out(index_real,3); % 交点的Z坐标（低点）
                n_receive=sum((z_jiao2>=80-4 & z_jiao2<=80+4));
                n_jie=n_jie+n_receive; % 总共有多少光线被接收
            end
        end
        eff_trun(i)=n_jie/((nsigma*ntau+1)*((dx+1)*(dy+1)-num_shadow)); % 计算i镜截断效率
    end
    %% 计算光学效率与热输出功率
    eff=eff_at.*eff_cos*eff_ref.*eff_sb.*eff_trun;%光学效率
    mean_eff=mean(eff);
    mean_eff_cos=mean(eff_cos);
    mean_eff_sb=mean(eff_sb);
    mean_eff_trunc=mean(eff_trun);
    eff_t(:,temp_t)=eff;
    power_t(temp_t)=dni*mir_wide*mir_length*sum(eff);%热功率
    dni_t(temp_t)=dni;
end
mean_dni=mean(dni_t);
eff_mean=mean(eff_t,2); % 计算每个镜一天的平均光学效率
[sort_eff_mean,index]=sort(eff_mean); % 对每个镜子平均光学效率进行排序，由低到高，下面要删除低效率的
for i=1:length(index)
    new_power=power_t-dni_t*mir_wide*mir_length*sum(sort_eff_mean(1:i)); % 把前ii个低效率的镜子去掉，功率所以要减去,如果全部去掉，功率变成0
    if mean(new_power)<3e4  %只算了一半，所以是3e4kW也就是30MW
        index_delete=i-1;   % 此时找到了要删除的最后一片镜子的序号！
        break;
    end
end

%% 绘图
id=1:length(X);
ind_new=setdiff(id,index(1:index_delete)); % 只要展示没有删除的点
figure
plot(X(ind_new),Y(ind_new),'o')
axis equal
figure
scatter(X(ind_new),Y(ind_new),20*eff_mean(ind_new),eff_mean(ind_new));
hold on
scatter(-X(ind_new),Y(ind_new),20*eff_mean(ind_new),eff_mean(ind_new));
axis equal
colormap jet
colorbar
xlabel('x')
ylabel('y')
title('最优化定日场布局以及平均光学效率分布图')

end
